#try to read in exonerate output and parse into gff3 format


#work on 1405-1 example, note no gaps in this

setwd('/Paterson/Datafiles/grouse')

exonerate.gff3 <- function(exonerate.file,fasta.file,gene_name,out.file='test_files/test.gff3',tmp.file="test_files/tmp.txt"){
	ex_out <- readLines(exonerate.file)
	
	writeLines("##gff-version 3",con=out.file)
	
	gff_out <- data.frame(start.line=grep('gff-version 2',ex_out),stop.line=grep('END OF GFF DUMP',ex_out))
	
	#start loop for every gene in gff_out
	for(i in 1:dim(gff_out)[1]){
		
		writeLines(ex_out[(gff_out$start.line[i]+8):(gff_out$stop.line[i]-2)],con=tmp.file)
	
		gff_gene <- read.table(tmp.file,col.names=c("seqname", "source", "feature", "start", "end", "score", "strand", "frame", "attributes"),sep="\t",stringsAsFactors=F,header=F)
		gff_gene2 <- gff_gene

		gff_gene2$seqname <- gene_name
		gff_gene2$feature[gff_gene$feature=="cds"] <- "CDS"
		gff_gene2$feature[gff_gene$feature=="splice5"] <- "five_prime_cis_splice_site"
		gff_gene2$feature[gff_gene$feature=="splice3"] <- "three_prime_cis_splice_site"

		gff_gene2$frame[gff_gene$feature=="cds"&gff_gene$frame=="."] <- "0"

		#annotate gene
		tmp <- unlist(strsplit(gff_gene$attributes[gff_gene$feature=="gene"][1]," ; "))
		if(i <10){
			tmp.parentID <- paste("gene0000",as.character(i),sep="")
			}else{
			tmp.parentID <- paste("gene000",as.character(i),sep="")
			}

		gff_gene2$attributes[gff_gene$feature=="gene"][1] <- 
		paste("ID=",tmp.parentID,";Name=",gff_gene2$seqname[1],".gene",i,sep="") 
		tmp.parentID <- sub('gene','mRNA',tmp.parentID)
		#annotate cds
		tmp.grep <- grep('cds',gff_gene$feature)
		for(j in 1:length(tmp.grep)){
			if(j < 10){
				tmp.cdsnm <- paste("cds0000",j,sep="")
				}else{
				tmp.cdsnm <- paste("cds000",j,sep="")
				}
			gff_gene2$attributes[tmp.grep[j]] <- paste("ID=",tmp.cdsnm,";Parent=",tmp.parentID,";Name=",gff_gene2$seqname[1],".cds",i,"_",j,sep="")
			}

		#annotate exons
		tmp.grep <- grep('exon',gff_gene$feature)
		for(j in 1:length(tmp.grep)){
			if(j < 10){
				tmp.cdsnm <- paste("exon0000",j,sep="")
				}else{
				tmp.cdsnm <- paste("exon000",j,sep="")
				}
			gff_gene2$attributes[tmp.grep[j]] <- 
			paste("ID=",tmp.cdsnm,";Parent=",tmp.parentID,";Name=",gff_gene2$seqname[1],".exon",i,"_",j,sep="")
			}

		#annotate introns
		tmp.grep <- grep('intron',gff_gene$feature)
		for(j in 1:length(tmp.grep)){
			if(j < 10){
				tmp.intnm <- paste("intron0000",j,sep="")
				}else{
				tmp.intnm <- paste("intron000",j,sep="")
				}
			gff_gene2$attributes[tmp.grep[j]] <- 
				paste("ID=",tmp.intnm,";Parent=",tmp.parentID,";Name=",gff_gene2$seqname[1],".intron",i,"_",j,sep="")
			#annotate 5' ss
			if(gff_gene$feature[tmp.grep[j]-1]=="splice5"){
				if(j < 10){
					tmp.ssnm <- paste("splicefive0000",j,sep="")
					}else{
					tmp.ssnm <- paste("splicefive000",j,sep="")
					}
		
				gff_gene2$attributes[tmp.grep[j]-1] <- 
					paste("ID=",tmp.intnm,";Parent=",tmp.intnm,";Name=",gff_gene2$seqname[1],".5splice",i,"_",j,sep="")
				}
	
			#annotate 3' ss
			if(gff_gene$feature[tmp.grep[j]+1]=="splice3"){
				if(j < 10){
					tmp.ssnm <- paste("splicethree0000",j,sep="")
					}else{
					tmp.ssnm <- paste("splicethree000",j,sep="")
					}
		
				gff_gene2$attributes[tmp.grep[j]+1] <- 
					paste("ID=",tmp.intnm,";Parent=",tmp.intnm,";Name=",gff_gene2$seqname[1],".3splice",i,"_",j,sep="")
				}

	
			}
		gff_gene3 <- gff_gene2[1,]
		gff_gene3[2,] <- gff_gene2[1,]
		gff_gene3[2,'feature'] <- 'mRNA'
		gff_gene3[2,'attributes'] <- gsub('gene','mRNA',gff_gene3[1,'attributes'])
		gff_gene3[2,'attributes'] <- paste(gff_gene3[2,'attributes'],";Parent=",sub('mRNA','gene',tmp.parentID),sep="")
		gff_gene3 <- rbind(gff_gene3,gff_gene2[-1,])

	
		#writeLines(paste("sequence-region")) #finish this later
		write.table(gff_gene3,file=out.file,append=T,sep="\t",col.names=F,row.names=F,quote=F)

		
		
		}
	cat("##FASTA\n",file=out.file,append=T)
	tmp.s <- strReadFasta(file=fasta.file)
	strWriteFasta(s=tmp.s,file=out.file,append=T,name=gene_name)
	cat("annotation written to ",out.file,"\n")
	
	}

#exonerate.gff3('test_files/test_prot2genome.out',fasta.file='test_files/cDNA_1405-1_matched.fasta',gene_name="cDNA_1405-1",out.file="test_files/test_fn.gff3")
#apollo complains
#sed '/intron/d' test.gff3 > test2.gff3

#try something similar to above but merge exons
get.exonerate <- function(exonerate.file,gene_name,tmp.file="test_files/tmp.txt"){
	ex_out <- readLines(exonerate.file)
	
	
	gff_out <- data.frame(start.line=grep('gff-version 2',ex_out),stop.line=grep('END OF GFF DUMP',ex_out))
	if(dim(gff_out)[1]==0){
		cat('no gff data found for ',exonerate.file,'\n')
		return(0)
		}
	#start loop for every gene in gff_out
	gff_gene2 <- list()
	for(i in 1:dim(gff_out)[1]){
		
		writeLines(ex_out[(gff_out$start.line[i]+8):(gff_out$stop.line[i]-2)],con=tmp.file)
	
		gff_gene <- read.table(tmp.file,col.names=c("seqname", "source", "feature", "start", "end", "score", "strand", "frame", "attributes"),sep="\t",stringsAsFactors=F,header=F)
		gff_gene2[[i]] <- gff_gene

		gff_gene2[[i]]$seqname <- gene_name
		gff_gene2[[i]]$feature[gff_gene$feature=="cds"] <- "CDS"
		gff_gene2[[i]]$feature[gff_gene$feature=="splice5"] <- "five_prime_cis_splice_site"
		gff_gene2[[i]]$feature[gff_gene$feature=="splice3"] <- "three_prime_cis_splice_site"

		gff_gene2[[i]]$frame[gff_gene$feature=="cds"&gff_gene$frame=="."] <- "0"

		#annotate gene
		tmp <- unlist(strsplit(gff_gene$attributes[gff_gene$feature=="gene"][1]," ; "))
		if(i <10){
			tmp.parentID <- paste("gene0000",as.character(i),sep="")
			}else{
			tmp.parentID <- paste("gene000",as.character(i),sep="")
			}

		gff_gene2[[i]]$attributes[gff_gene$feature=="gene"][1] <- 
		paste("ID=",tmp.parentID,";Name=",gff_gene2[[i]]$seqname[1],".gene",i,sep="") 
		tmp.parentID <- sub('gene','mRNA',tmp.parentID)
		#annotate cds
		tmp.grep <- grep('cds',gff_gene$feature)
		for(j in 1:length(tmp.grep)){
			if(j < 10){
				tmp.cdsnm <- paste("cds0000",i,sep="")
				}else{
				tmp.cdsnm <- paste("cds000",i,sep="")
				}
			gff_gene2[[i]]$attributes[tmp.grep[j]] <- paste("ID=",tmp.cdsnm,";Parent=",tmp.parentID,";Name=",gff_gene2[[i]]$seqname[1],".cds",i,"_",j,sep="")
			}

		#annotate exons
		tmp.grep <- grep('exon',gff_gene$feature)
		if(length(tmp.grep)>0){
			for(j in 1:length(tmp.grep)){
				if(j < 10){
					tmp.cdsnm <- paste("exon0000",j,sep="")
					}else{
					tmp.cdsnm <- paste("exon000",j,sep="")
					}
				gff_gene2[[i]]$attributes[tmp.grep[j]] <- 
				paste("ID=",tmp.cdsnm,";Parent=",tmp.parentID,";Name=",gff_gene2[[i]]$seqname[1],".exon",i,"_",j,sep="")
				}

			}
		
		#annotate introns
		tmp.grep <- grep('intron',gff_gene$feature)
		if(length(tmp.grep)>0){
			for(j in 1:length(tmp.grep)){
				if(j < 10){
					tmp.intnm <- paste("intron0000",j,sep="")
					}else{
					tmp.intnm <- paste("intron000",j,sep="")
					}
				gff_gene2[[i]]$attributes[tmp.grep[j]] <- 
					paste("ID=",tmp.intnm,";Parent=",tmp.parentID,";Name=",gff_gene2[[i]]$seqname[1],".intron",i,"_",j,sep="")
				#annotate 5' ss
				if(gff_gene$feature[tmp.grep[j]-1]=="splice5"){
					if(j < 10){
						tmp.ssnm <- paste("splicefive0000",j,sep="")
						}else{
						tmp.ssnm <- paste("splicefive000",j,sep="")
						}
		
					gff_gene2[[i]]$attributes[tmp.grep[j]-1] <- 
						paste("ID=",tmp.intnm,";Parent=",tmp.intnm,";Name=",gff_gene2[[i]]$seqname[1],".5splice",i,"_",j,sep="")
					}
	
				#annotate 3' ss
				if(gff_gene$feature[tmp.grep[j]+1]=="splice3"){
					if(j < 10){
						tmp.ssnm <- paste("splicethree0000",j,sep="")
						}else{
						tmp.ssnm <- paste("splicethree000",j,sep="")
						}
		
					gff_gene2[[i]]$attributes[tmp.grep[j]+1] <- 
						paste("ID=",tmp.intnm,";Parent=",tmp.intnm,";Name=",gff_gene2[[i]]$seqname[1],".3splice",i,"_",j,sep="")
					}

	
				}
			}

		}
	gff_gene2
	}
#debug(get.exonerate)
tst <- get.exonerate('test_files/test_prot2genome.out',gene_name='cDNA1405-1')

tst2 <- tst[[1]][grep('CDS',tst[[1]]$feature),]
for(i in 2:length(tst)){
	tst2 <- rbind(tst2,tst[[i]][grep('CDS',tst[[i]]$feature),])
	}

add.CDS <- function(embl.file,cds.list,out.file,sequence=NULL){
	#insert cds info into an embl file
	#still haven't worked out what to do with gaps
	
	embl.lines <- readLines(embl.file)
	embl.seq <- strReadEmbl(embl.file)
	sq.start <- grep('^SQ',embl.lines)
	writeLines(embl.lines[1:(sq.start-2)],out.file)
	for(i in 1:length(cds.list)){
		tmp.cds <- cds.list[[i]][grep('CDS',cds.list[[i]]$feature),]
		writeEmblLine(out.file,code="FT",header="CDS",
			text=paste("join(",paste(paste(tmp.cds$start,'..',tmp.cds$end,sep=""),collapse=","),")",sep=""),
			nextfield=FALSE)
		tmp.cdna <- character()
		for(j in 1:dim(tmp.cds)[1]) tmp.cdna <- paste(tmp.cdna,substr(embl.seq,tmp.cds$start[j],tmp.cds$end[j]))
		tmp.cdna <- gsub(" ","",tmp.cdna)
		writeEmblLine(out.file,code="FT",text=paste('/product=\"',paste(tmp.cds$seqname[1],"cds",i,sep="_"),'\"',sep=""),nextfield=FALSE)
		writeEmblLine(out.file,code="FT",text=paste('/note=\"',paste("score", cds.list[[i]]$score[1],sep=" "),'\"',sep=""),nextfield=FALSE)
		writeEmblLine(out.file,code="FT",text=paste('/translation=\"',
			strTranslate(tmp.cdna),'\"',sep=""),nextfield=FALSE)
				}
	tmp.con <- file(out.file,open="a")
	writeLines(embl.lines[(sq.start-1):length(embl.lines)],tmp.con)
	close(tmp.con)
	cat("annotation written to ",out.file,"\n")
	
	}	
#add.CDS("test_files/cDNA_1405-1_matched.embl",cds.list=tst,out.file="test_files/test2.embl")

